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Collisionless evaporation from cluster elliptical galaxies: 
a contributor to the intracluster stellar population 
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Abstract. By means of simple numerical models we discuss whether collisionless stellar evaporation from cluster 
elliptical galaxies could be an effective mechanism for the production of intracluster stellar populations. The 
effectiveness of this mechanism is due to the fact that, for realistic galaxy and cluster models, the galaxy oscillation 
periods near equilibrium configurations in the cluster tidal field are of the same order of stellar orbital times in 
the external parts of the galaxies themselves. With the aid of Monte-Carlo simulations we explore the evolution 
of stellar orbits in oscillating galaxies placed near different equilibrium positions. We found that, over an Hubble 
time, the main effect is a substantial expansion of the galactic outskirts, particularly affecting the galaxy at the 
cluster center and those orbiting near the cluster core radius: overall, approximately the 10% of the galaxy mass 
is affected. Thus, the proposed mechanism could be of some importance in the shaping of the halo of cD galaxies 
and in making easier "galaxy harassment" in the formation of the intracluster stellar population. 
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1. Introduction 

Observational evidences of an intracluster stellar popula- 
tion (hereafter ISP) are mainly based on the detection 
of planetary nebulae and red giant branch stars freely 
floating in the intergalactic space. For example, Theuns & 
Warren (1996) identified 10 intergalactic planetary nebu- 
lae in the Fornax cluster, while Mendez et al. (1997) de- 
tected 11 intergalactic objects in the Virgo cluster whose 
cumulative luminosity functions is in good agreement 
with planetary nebula luminosity function. In addition, 
Ferguson et al. (1998) identified an intergalactic red gi- 
ant branch stars population in the Virgo cluster, while 
Feldmeier et al. (1998) (with observations of three blank 
fields in the Virgo cluster), confirmed that a significant 
fraction of Virgo's starlight is due to the ISP. More re- 
cently Okamura et al. (2002) have identified 38 candidates 
of intracluster planetary nebulae in the core of the Virgo 
cluster. Overall, the data suggest that approximately 10% 
(or even more) of the stellar mass of the cluster is in inter- 
galactic stars (e.g., see Ferguson et al. 1998, Arnaboldi et 
al. 2002, Durrel et al. 2002, Arnaboldi et al. 2003, Totani 
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2003, see also Arnaboldi, Gerhard & Freeman 2003, and 
references therein). 

The usual scenario assumed to explain the finding 
above is that the tidal interactions between galaxies (for 
example during a fast encounter, e.g., see Merritt 1983, 
1984, 1985), and of galaxies with the cluster gravitational 
field, lead to a substantial stripping of stars from galaxies 
to the parent cluster (the so-called "galaxy harassment" 
scenario, see, e.g., Moore et al. 1996; Napolitano et al. 
2003). In this paper we explore an additional "stripping" 
mechanism, namely we quantitatively discuss the effect 
of resonant interaction between stellar orbits inside the 
galaxies and the cluster tidal field (hereafter, CTF). The 
present study is supported by the fact that the charac- 
teristic oscillation times of a galaxy near its equilibrium 
position in the CTF and the mean stellar orbital times 
in the galaxy external regions are of the same order of 
magnitude, as already recognized by Ciotti & Giampieri 
(1998, hereafter CG98). 

In fact, Hawley & Peebles (1975) reported a possible 
indication (confirmed by Thompson 1976), that the galax- 
ies are preferentially aligned along the radius vector to the 
center of the cluster for the Coma cluster; they also sug- 
gested the CTF as the possible cause of the alignment. The 
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best evidence for alignment is for brightest cluster galax- 
ies: Adams, Strom & Strom (1980) found, in 7 very elon- 
gated clusters, a general trend for ellipticals to be aligned 
with the cluster major axis, and also Trevese, Cirmele 
& Flin (1992) found a strong alignment of the bright- 
est galaxy major with the long axis of the parent clus- 
ter. Numerical N-body simulations (Ciotti & Dutta 1994, 
hereafter CD94) indeed confirmed the hypothesis that the 
CTF could be at the origin of the observed alignment and, 
in particular, revealed that their model elliptical galaxies 
(Es) behaved with good approximation as rigid bodies. In 
a complementary analytical exploration of this problem 
CG98 determined the equilibrium positions of triaxial el- 
lipsoids in the CTF for various cases, and showed that 
the oscillations period of the galaxies, when slightly dis- 
placed from their equilibrium configurations, are of the 
same order of magnitude of the stellar orbital periods in 
the galaxy outskirts. This curious finding naturally leads 
to ask what could be the effect of possible "resonances" 
between stellar orbital periods in cluster Es and their os- 
cillation periods. Unfortunately, the A-body simulations 
of CD94 were characterized by a limited number of parti- 
cles (i.e., N ~ 3 x 10 4 ) and by the usual softening, so that 
evaporations rates could not be properly investigated. 

In order to investigate the scenario described above, 
Muccione & Ciotti (2003a,b, hereafter MC03a,b) per- 
formed a few preliminary Monte-Carlo simulations of or- 
bital exploration of oscillating galaxies, and the results 
were encouraging, with estimated escape rates summing 
up to 10 per cent of the initial galaxy mass over an Hubble 
time. In the present paper we considerably extend these 
two previous investigations, by assuming that the galaxies 
arc triaxial ellipsoids and two cases are considered: when 
the center of mass of the ellipsoid is at rest at the equilib- 
rium point of the field generated by the cluster, and when 
it is placed on a circular orbit around the center of a spher- 
ically symmetric cluster. After deriving the equations of 
the motion of a star inside an oscillating galaxy, we in- 
tegrate numerically these equations under different initial 
conditions realized by using a Monte-Carlo extraction. 

The paper is organized as follows. In Sect. 2 we briefly 
review the proper physical setting of the problem, and in 
Sect. 3 we describe the adopted galaxy and cluster models. 
In Sect. 4 we describe in the numerical integration scheme, 
while in Sect. 5 we present the main results, that are finally 
summarized and discussed in Sect. 6. In the Appendix, 
the specific technique adopted to obtain the gravitational 
field inside the galaxy models and other useful dynamical 
quantities are briefly described. 



2. The physical setting 

2.1. Triaxial galaxy at the center of a triaxial cluster 

Following the analytical treatment of CG98, we start pre- 
senting the simple case of a spinless galaxy with its center 
of mass at rest at the center of a triaxial galaxy clus- 
ter. In CG98 it was shown that the galactic equilibrium 



configurations correspond to the galaxy inertia ellipsoid 
oriented along the CTF principal directions: without loss 
of generality we assume that in the (incrtial) Cartesian co- 
ordinate system C, the CTF tensor T is in diagonal form, 
with components T t (i = 1,2,3). 

By using three successive counterclockwise rotations 
(ip around x axis, d around y' axis and ip around z" axis), 
CG98 showed that the linearized equations of motion for 
the galaxy near the equilibrium configurations can be writ- 
ten as 



.. AT 32 A/ 32 

<P = r 



h 

AT 31 A/ 31 



0. 



h 

IP = t V', 



(1) 



where ATy = Ti—Tj, and 7j are the principal components 
of the galaxy inertia tensor. If we also assume that T\ > 
Ti > T 3 and I\ < 7 2 < I3, then the equilibrium position 
of eq. (1) is stable (for the geometrical interpretation of 
this condition see Sect. 3), and the solution of eq. (1) is 

ip = (fM cos(u) v t), -d — § M cos(uj#t), ip — ipM cos(o^i), (2) 

where 



AT 23 A/ 32 



h 



AT 13 A/ 31 



AT 12 A/ 2 



(3) 



are the three independent oscillations frequencies. For sim- 
plicity we assume that at t = the galaxy is at the max- 
imum deviation from the equilibrium configuration (with 
null angular velocity). 

For computational reasons the best reference system 
in which calculate and discuss the stellar orbits is the os- 
cillating reference system C in which the galaxy is at rest, 
and its inertia tensor is in diagonal form. As a consequence 
of this choice, the equation of the motion for a galactic star 
in C is 



I = TZ T x - 2ft A v - ft A £ - ft A (ft A £), 



(4) 



where £ and v = £ are the stellar position and veloc- 
ity in C", the suffix "T" means transpose, and 1Z is the 
orthogonal transformation matrix between C and C so 
that x = 1Z£. The explicit form of 1Z is given by eq. (9) 
in CG98: as well known, this matrix define the angular 
velocity of C that, expressed in C, reads 



ft = (v3 cos $ cos -0 + i?sun/>, 

—tp cosz9 sin ^ + § cos?/;, y>sini9 + ip). 



(5) 



If we indicate with g = g (£) the galactic gravitational 
potential in C, then in cq. (4) 

TZ T i = -Ve(f> g + K T T1l£, (6) 
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where is the gradient operator in C, and the direct 
effect of the cluster gravitational field is obtained in the 
tidal approximation, fn this scheme, each star experience 
two different influences from the ambient cluster: an in- 
direct one, due to the induced oscillatory motion of the 
galaxy, and the direct acceleration due to the galaxy and 
the cluster mass distributions. Of these two last terms, 
the by far more important within the galaxy is the galac- 
tic gravitational field, fn Sects. 3 and 4, and in Appendices 
B and C, we give the explicit expression for the quantities 
appearing in eqs. (4) and (6) for the adopted galaxy and 
cluster models. 



2.2. Triaxial galaxy on circular orbit in a spherical 
cluster 

The second, more complicate (but still highly idealized, 
although more astrophysically relevant) case of interest is 
represented by a triaxial galaxy with its center of mass 
in circular orbit of radius r in a plane (say, z = 0) of a 
spherical cluster, with angular velocity f2 c = (0,0, fi c ), 
where 



GM c (r) 



(7) 



and where M c (r) is the cluster mass within r. Following 
the treatment of CG98, the analogues of eqs. (1) are 



<p = — 



n 2 c Ai 32 , n c (h - ai 32 ) . 



h 



-<p- 



h 



(O c 2 + AT 13 )A7 3 i ^ n c (/ 2 - A/ 3 i) 



(8) 



while the equation for ip remains unchanged. The explicit 
solution of eqs. (8) is given in Appendix A. fn this case 
there are two different equilibrium configurations: the first 
corresponds to the galaxy's major axis directed toward 
the cluster center, and the second to the galaxy major 
axis perpendicular to the orbital plane (see CG98). For 
simplicity in this paper we restrict to the first case only. 

Due to the system configuration we now need two co- 
ordinate transformations in order to obtain the equations 
of motion for a star in the galaxy. The first transforms 
the (cluster) inertial system C in the system C c centered 
on the galaxy center of mass (rotating with the angular 
velocity fl c and with its z axis parallel to the z axis of C) , 
and the second transforms C c in the system C in which 
the galaxy is at rest and its inertia tensor is in diagonal 
form. The relation between the position vector x in C and 
the position vector £ in C is given byx = r + 7£'£, where 
1Z' = 1Z C 1Z. 1Z is the same as in eq. (4), while 



' cosCl c t — sinfi c t 
1Z C = | sinf2 c i cos fl c t 
f 



(9) 



Straightforward algebra then shows that the equations of 
motion for a star as seen from C are similar to eq. (4), 



where f2 is substituted by f2' — ft + TZ T fl c , f2 is given in 
eq. (5), and 

lZ T n c = (sin tp simp — cos </? sin $ cost/*, 

sin ip cos tp + cos ip sin t? sin tp, cos ip cos $)fi c .(10) 

Moreover, 7\L T x is substituted by 7£' T (x — r) ~ — V^^> g + 
1Z T T1Z£, where now T is the CTF tensor as seen in C c . 
Note that in tidal approximation, the center of mass of 
the galaxy will consistently remains in the circular orbit, 
independently of the galaxy oscillations (see CG98). In 
Sects. 3 and 4, and in Appendices A, B, and C we give 
the explicit expression for the various quantities of interest 
for the adopted galaxy and cluster models. 



3. Galaxy and cluster models 

We now present the specific galaxy and cluster density 
distributions adopted in the numerical simulations, and 
we derive the required explicit expressions for the galaxy 
inertia tensor and the CTF: a similar approach can be 
found in Valluri (1993). For simplicity we assume that the 
galaxy and cluster densities are stratified on homeoids. In 
MC03ab we adopted triaxial n = Ferrers (1887) and 
Hcrnquist (1990) models, respectively: here we recall the 
results obtained when adopting the ellipsoidal generaliza- 
tion of the widely used 7- models (Dehnen 1993): 



r(m) 



M g 3- 7 



1 



aia 2 «3 47r mT(l + to) 4 ~t ' 



0< 7 <3, (11) 



while we perform the extended exploration of the param- 
eter space by using the density distribution 



p g (m) = 



M s 15 



a.\ctict 3 2ir (1 + m) 6 ' 
where M g is the total mass of the galaxy, and 



3 

E 

i=l 



£ 2 



a 



(12) 



(13) 



Quite obviously, the density profiles in eq. (11) are more 
realistic than that in eq. (12): however, for technical rea- 
sons described in Appendix C, this latter density profile 
allows for much more faster numerical simulations (i.e., a 
larger number of particles). Fortunately, the obtained es- 
cape rates in all the explored models, here and in MC03ab, 
are remarkably similar, and so can be considered quite ro- 
bust (within the considered scenario). 

The inertia tensor of generic homeoids arc given by 



7- 47T 2 2 



(14) 



where h g = J °° p g (m)m 4 dm, and so the stability require- 
ment Ji < I2 < I3 is satisfied. Note also that, according 
to eq. (3), the frequencies for homeoidal stratifications do 
not depend on the specific density distribution assumed, 
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but only on the quantities (ai, 02,0:3). For a better com- 
parison with observations, in the following we use the two 
cllipticities 1 



— = 1 — e, — = 1 — 77, e < 77 < 1. 

Oil Ot\ 



(15) 



A rough estimate of characteristic stellar orbital times 
inside m for an homeoidal distribution is given by 
-Porb(^) — 4P dyn (m) = ^37r/Gp g (m), where for a 7 
model the mean galaxy density inside m is 



p g (m) 



3M e 



1 



47raia 2 a 3 mT(l + to) 3 ~t ' 
and the galaxy mass inside m is given by: 



(16) 



f m m 3 ~ 7 
M g (m) = 4TOia 2 a 3 y p g (t)t 2 (it = M s - _ 7 . (17) 



Thus, 



P orb (m) ~ 9.35 x 10' 




Ul-eXl-r?) 



m 7/2 (l + m) (3 ^ )/2 yrs, 



(18) 



where M gi n = M g /10 11 M Q and a^i = ai/kpc is the 
galaxy "core" major axis 2 ; thus, in the outskirts of normal 
galaxies orbital times well exceeds 10 8 or even 10 9 yrs. A 
similar conclusion is reached also for the model in eq. (12). 

To describe the cluster density distribution we use the 
simple mass profile 



Pc(m) 



PcO 



(1 + m 2 ) 2 ' 



(19) 



where m is given by an identity similar to eq. (13), with 
ai > a 2 > 03; in the case of the spherical cluster we 
assume a\ —02 = 03. In the two following Sections we 
will compute the CTF components associated with eq. 
(19). 



p < v < 1, and the expansion of the quantities Wi for 
small flattenings is given in Appendix C. 

In order to determine the galaxy oscillation frequencies 
we need, according to eqs. (3) and (20), to have a realis- 
tic estimate of p c0 . Unfortunately this quantity is not well 
measured in real clusters, and for its determination we use 
the virial theorem, M c a\ — —U (U is the cluster gravi- 
tational potential energy and <7y is the virial velocity dis- 
persion, that we assume to be estimated by the observed 
velocity dispersion of galaxies in the cluster). The explicit 
calculations of the factor 27rG / o c o is presented in Appendix 
B. 

We now compare the galactic oscillation periods P v = 
2tt/ui ip , P$ — 27r/a>,j, and P^ = with the charac- 

teristic stellar orbital times P rb in galaxies: for simplicity 
we give here the expansions for small cluster and galac- 
tic flattenings. From eqs. (3), (14), (20), and using eqs. 
(C.8)-(C.12) wc finally obtains to the leading order in the 
flattenings 



8.58 x 10 8 



fll.250 



\J{v — p)(v - e) OV4000 



yrs, 



P<> * 



Pip — 



8.58 x 10 s ai 251 



vil °V,iooo 
8.58 x 10 8 oi 25 



/pie 



250 

fv.iooo 



yrs, 



yrs, 



(21) 
(22) 
(23) 



where we normalized cr v to 1000 km/s and ai to 250 
kpc. So, a comparison of quantities above with eq. (18) 
shows that in the outer halo of giant Es the stellar orbital 
times can be of the same order of magnitude as the oscil- 
latory periods of the galaxies themselves. For example, in 
a relatively small (Hcrnquist) galaxy of M gi n =0.1 and 
ais — 1, P rb — 1 Gyr at m ~ 10 (i.e., at ~ 5R e ), while 
for a galaxy with M gi n = 1 and ai,i = 3 the same orbital 
time is at m ~ 7 (i.e., at ~ 3.5i? c ). 



3.1. Galaxy at the cluster center 

In Appendix B it is shown that the tidal field compo- 
nents at the center of a non-singular homeoidal distribu- 
tion (such as that in eq. [19]), are given by 



T 4 = -2nGp C QW l . 



(20) 



Note that the quantities w; do not depend on the specific 
density profile, and that w\ < W2 < W3 for oi > a 2 > a 3 , 
thus fulfilling the conditions for stable equilibrium in eq. 
(1). In analogy with definitions (15), we introduce the two 
cluster cllipticities, a 2 /ai = 1 — p and 03/01 = 1 — u, with 



1 In fact, it can be easily proved the for generic homeoidal 
distributions, e and n defined according to eq. (15) are the true 
ellipticities in the relative projection planes. 

2 We recall here that for the spherically symmetric Hernquist 
model, (i.e., Qi = 02 = 03 in eq. [12]), _R C ~ 1.8ai, while for 
the density distribution in eq. (12) R c ~ 0.75qi. 



3.2. Galaxy in circular orbit 

In the case of spherical cluster, the CTF tensor in the 
reference system C c introduced in Sect. 2.2 is 



'3q- 2 s 
T = -n^(r) I 1 I , 
1 



(24) 



where q(r) = p c (r)/p c (r) and p c = 3M c (r)/47rr 3 (see, e.g., 
CD94, CG98). From eq. (19) wc obtain 

M c (r) = 27ra?/9 c o ^arctanf - - ^-^ j , (25) 

(26) 



ft( r ) = ^F (arctanf- 



. . 2f 3 
q(r) = I arctanr 



3(l + f 2 ) 2 



l + f 



P2 



(27) 
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where f = r/a\. In analogy with Sect. 3.1, we now com- 
pare the galaxy oscillation periods with the mean stellar 
orbital times, and for simplicity we restrict to the special 
case of the galaxy oscillating around the ^3-axis only, i.e., 
we fix ip = and # = in eq. (10). Thus, from the results 
of Sect. 2.2, and from eq. (24) 



p 27T 



2tt 



1 



wv- ^c(r) y/3e[l - q(r)} 



(28) 



where we used the last of eqs. (3) and the expansion for 
small flattcnings of the coefficient involving the galaxy mo- 
ments of inertia. In Fig. 1 we plot P$ as a function of r 
for the two spherical cluster models used in the simula- 
tions described in Sect. 5.2, in the cases of an E2 and and 
E4 galaxy; note that P^ decreases at increasing galaxy 
flattening, while its radial trend is non monotonic, with 
a minimum at r ~ a\, the cluster core radius. The pecu- 
liarity of this location has been already noted (e.g., see 
CD94), and its impact on orbital evolution will be evident 
in Sect. 5.2. 




Fig. 1. The oscillation period around the £3 || fl c axis for a 
galaxy in circular orbit of radius r in a cluster of core radius 
ai and velocity dispersion av- Solid lines refer to an E2 galaxy, 
while dotted lines to an E4 galaxy. Cluster model Ca (heavy 
lines) is characterized by a\ — 250 kpc and av = 1000 km/s, 
while model Cb by a\ — 100 kpc and av = 500 km/s. 



4. The numerical integration scheme 

4.1. Dimensionless equations 

In the numerical simulations we use as time normalization 
the Hubble time tu = 15 Gyrs, as mass normalization 



the galaxy mass M g , and as length normalization ot\, the 
homeoidal scale-length corresponding to the galaxy major 
axis. Accordingly eq. (4) in dimensionless form reads 



dr 2 



-2uj c t n Sl Av — cj^[SI A £ + ft A (0 A £)], (29) 



note that from 
and so from eqs. 



where r = t/% & n d g = {GM g /a\ 
eq. (3) u>i = uj c &i where to, 
(2) and (5) SI = u) c Sl. In addition, consistently with the 
adopted "small oscillations approximation" , the angular 
velocity is used in linearized form 



(30) 



i.e., in eq. (5) we retained only linear terms in tfM, $m> 
V'm! note that we neglect the direct effect of the cluster 
gravitational field on stellar orbits (cf. with eq. [6]). 

According to the discussion in Sect. 2.2, the dimen- 
sionless equations of the motion for a star in the galaxy 
on a circular orbit arc similar to eq. (29), where now 
lo c = y/2irGp c , we again disregard the direct effect of the 
CTF, SI is substituted by SI' = SI' /u> c , and to the linear 
order 

n'~ {<p-m c ,# + ipn c ,i> + n c ). (31) 

Note that Sl c /u; c = yplji. 
4.2. Initial conditions 

To integrate the second order differential equations (29) 
describing the motion of a test star in the non-inertial ref- 
erence frame fixed in the galaxy, we use a code based on an 
adapted f77, double-precision 4th order Rungc-Kutta rou- 
tine (Press et al. 1986), with adaptive time step. In prac- 
tice, for a given galaxy and cluster model, we analyzed the 
behaviour of a large number of orbits (Atot ~ 10 4 -j- 10 6 ) 
over an equivalent Hubble time: from this point of view 
our simulations are nothing else that a large collection 
of one-body problems in time dependent potentials. Of 
course, the initial conditions for each orbit must be re- 
alized with some care. First of all, in order to properly 
sample the galaxy density profile the initial positions are 
obtained from the Von Neumman rejection method (Press 
ct al. 1986) which uses the (elliptical) cumulative mass 
distribution M g (m)/M g as probability function for the 
homeoidal radius; the two position angles are then ex- 
tracted from uniformly populated standard angular inter- 
vals (i.e., < ip < 2tt and — 1 < coa-d < 1). 

A more delicate issue is represented by the assignment 
of the three components of the initial velocity. In principle, 
one should perform a von Neumann rejection method in 
phase-space, provided the model phase-space distribution 
function is known at each point. Unfortunately, for tri- 
axial galaxies the only sufficiently general method is the 
Schwarzschild (1979) orbit linear superposition method, 
which however for our purposes would be very demanding 
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from the numerical point of view. Here we adopt a more 
empirical approach, whose consistency is however readily 
verified. In some simulations, for each initial position in 
configuration space, we calculate the 1-dimensional veloc- 
ity dispersion a 1 from the Jeans equations of the spherical 
analogue of the density distribution in cqs. (11) and (12), 
we choose a dimensionless parameter \ < 1; an d we ran- 
domly distribute over the three directions the initial par- 
ticle velocity so that its square modulus equals 3xcr 2 . In 
other simulations, we same procedure was repeated, but 
assuming a fixed (within a given simulation) fraction \ of 
the local escape velocity v esc . 

For each simulation (i.e., for assigned galaxy and clus- 
ter models, maximum oscillation angles, and \ value) we 
actually run the code twice: in the first simulation the 
galaxy maximum oscillation angles are set to zero (in other 
words the galaxy is at rest at the equilibrium position in 
the CTF), and all the initial conditions leading to unac- 
ceptable orbits for the unperturbed case are discarded. 
In this way we can also test the code for energy con- 
servation. In all the computed models the "unperturbed" 
case pass very well the acceptability test, with relative er- 
rors on energy conservation for each orbit never exceeding 
~ 10~ 5 over an Hubble time, and without significant spu- 
rious evaporation. Thus, we can safely assume that the 
evaporation rates found in the explored models are a gen- 
uine product of galaxy oscillations, and no a product of 
and inconsistent assignment of the initial conditions. 

All the simulations have been performed on 
GRAVITOR, the Geneva Observatory 132 proces- 
sors Beowulf cluster 3 : for reference, the computation of 
10 4 orbits requires ~ 2 hours when using 10 nodes. 

5. Results 

5.1. Galaxy at the center of a triaxial cluster 

Here we describe the results of the first of the two cases 
discussed in Sect. 2, namely the case of a triaxial galaxy 
at the center of a triaxial cluster. In this case we already 
have preliminary informations for galaxy models described 
by Ferrers (1887, see MC03a) and by the more realistic 
Hernquist (see MC03b) density distribution. In this Sect, 
we complement these informations with the results ob- 
tained when using the density profile in eq. (12). However 
in the present for Hernquist models and at vari- 

ance with Ferrers distribution, the concept of "escape" is 
quite subtle. In fact, while in truncated galaxy models the 
escape criterion is obvious, in untruncated models the "es- 
caping" particles must be identified in a different way. In 
particular we choose to plot, for each test particle in a 
given model, the ratio m max /m; vs. m;, where m; is the 
ellipsoidal radius of a given initial condition, while m max 
is the maximum ellipsoidal distance reached by that parti- 
cle over an equivalent Hubble time. We present here only 
a representative simulation among the several performed, 

3 For technical detail on GRAVITOR see 
http://obswww.unige.ch/ pfennige/gravitor/gravitor _e.html 



because they remarkably confirm the previous results for 
Hernquist models described in MC03b. In particular, the 
shown simulation refers to a galaxy model characterized 
by the two flattenings e = 2/5 and r/ = 16/25; the same 
ellipticities have been also adopted for the cluster, i.e. 
/i = 2/5 and v = 16/25. Thus, the adopted galaxy model 
would be classified, as a function of the observation an- 
gle, from an EA up to an E6 galaxy. Its mass is fixed to 
M g = 1O 11 M , while its core semimajor axis is fixed to 
ct\ = 6.58 kpc (corresponding to an effective radius in the 
spherical model of ~ 5 kpc). The cluster is characterized 
by try = 1000 km/s and a\ = 250 kpc. Initial conditions 
were generated with mi < 8 and with a square modulus 
of the initial velocity of each particle equals to the local 
0.3v^ sc . The galaxy maximum oscillation angles along the 
three directions are lpm = $m = V'm = 0.25. 

The results are shown in Fig. 2. Note how in the unper- 
turbed case, with the adopted x value, ra max /TOi increases 
only slightly, while in the oscillating case this number 
reaches higher values, of the order of 1.5 or more. These 
results are very similar to those described in MC03b (see 
Figs. 1 and 2 there). Note also how the expansion effect 
is more important in the external galaxy regions (affect- 
ing particles with mj>3), even though we recall that in 
the simulations the direct effect of the CTF is not consid- 
ered: as a consequence, the obtained galaxy expansion is 
a genuine effect of galactic oscillations. By the way, this 
is a quite interesting result, that one could relate to the 
process of cD formation. In fact, it is well known that 
it is not easy to grow with galactic cannibalism, at least 
in numerical simulations, the extended cD halos (e.g., see 
Nipoti et al. 2003, and references therein). However, one 
should recall that the CTF in shallow cores is compressive 
(e.g., see Valluri 1993, CD94): as a consequence, the final 
galaxy density profile will be determined by the competing 
effect of galaxy oscillations and the cluster gravitational 
field. Thus, in order to quantitatively discuss the possi- 
bility that the galaxy sloshing at the cluster center is at 
the origin of the cD halos, high-resolution, self-consistent 
N-body simulations are required. 

5.2. Galaxy on circular orbit in a spherical cluster 

We now move to describe the results of the simulations of 
galaxies in circular orbit around the cluster center: altough 
still highly idealized, this cases apply to a larger number of 
galaxies than the previous one. At variance with the case 
of a galaxy at the cluster center, galaxies in circular orbit 
were not studied before, so here we present a larger set 
of simulations, organized in two groups. In the first set of 
simulations the spherical cluster is characterized by ay = 
1000 km s _1 and a\ — 250 kpc (model Ca), while in the 
second set ay = 500 km s _1 and a\ = 100 kpc (model Cb). 
For each model cluster we study three different positions 
of the galaxy center of mass, namely circular orbits with 
r/cii = 0.5 (models Cai and Cbi), with rja\ = 1 (models 
Cac and Cbc), and finally with r/ai = 2 (models Cao 
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Fig. 2. The ratio m max /mi vs. mi for the subset of particles 
with mi > 2, evaluated over tu- Upper panel: unperturbed 
galaxy with ipM = $m = ipM = 0. Lower panel: the same 
model with ipyi — $m = i>M = 0.25 rad. 



and Cbo). In order to reduce the dimensionality of the 
parameter space, in all the simulations we fix a galaxy 
model with M g = 10 n M Q , a x = 7 kpc, e = 0.2 and 
n = 0.4 (thus corresponding to a less flattened galaxy 
than that used in Sect. 5.1). Finally, we explore different 
cases of initial velocities of the test particles, with \ m 
the range 0.1 - 0.5. In all cases, we restrict our analysis 
to the case of a galaxy oscillating only around the £3 axis 
(the axis perpendicular to the orbital plane of the galaxy 
center of mass), and we adopt tpu = 0.2 rad. 

We note here that the set-up of the initial conditions is 
considerably more time-expensive that in the case of the 
galaxy at the cluster center. In fact, in the present case, 
"galaxy at rest" means that the galaxy center of mass 
is actually (uniformly) rotating around the cluster center, 
and this induces a centrifugal acceleration on the galaxy 
stars, as apparent by setting ip = i3 = ip = in cq. (31). 
Thus, if one erroneously adopt as stellar orbit initial condi- 
tions the same conditions adopted in the previous Section, 
the result is a vigorous stellar escape, just due to the cen- 
trifugal acceleration. Accordingly, test numerical simula- 
tions (not shown here) were characterized by a significant 
number of stars escaped up to m max ~ 100. Of course 
this is not a physical process, but only the result of the 
non-equilibrium initial conditions: in order to avoid this 
problem we proceeded as follows. For each galaxy model, 
we first performed the orbital calculation for a large num- 
ber (say 10 5 ) initial conditions arranged as in Sect. 5.1, 
maintaining the galaxy in uniform rotation around the 
cluster center and without oscillations. We then discarded 



all the initial conditions corresponding to orbits for which, 
over an Hubble time, m max /mi>1.2, in order to mimic the 
result shown in the bottom panel of Fig. 2. We finally used 
the remaining initial conditions to study the orbital evo- 
lution in the corresponding oscillating galaxy. The results 
are summarized in Fig. 3. 

Several interesting comments can be made from in- 
spection of Fig. 3. The first, and more obvious, is that 
also in the rotating galaxies, oscillations are expected to 
produce stellar evaporation. However, in the rotating cases 
we invariably found that the largest m max /TO; are larger 
that the corresponding quantity relative to the galaxies 
at the cluster center (cfr. Fig. 2), and this is due to the 
additional effect of the centrifugal term on escaping stars. 
Strictly related to this point is also the systematic be- 
havior of TOmax/wi shown in Fig. 3, where it is apparent 
how the largest ratios systematically decrease for galaxy 
models placed at r = a\, a\j1 and 2oi, respectively. This 
non-monotonic trend si just a reflection of the same be- 
havior of P^,, shown in Fig. 1: shorter the galaxy oscilla- 
tion times, larger are the TOmax/™. that can be reached by 
escaping stars. In this respect, the peculiar dynamical sit- 
uation of galaxies orbiting near the cluster core radius it 
is well known: not only there is a change in the topological 
nature of the CTF there (see, e.g., CD94), but also oscil- 
lation periods are the shortest (for a given galaxy model). 
The third comment about Fig. 3 is the common pattern 
of the radial behavior of the upper envelope of m max /mi, 
i.e., a quite well defined, sharp increase followed by a some- 
what smoother decline at increasing m\. This behavior ba- 
sically corresponds to the place in the galaxy where orbital 
times (see eq. [18]) become comparable with the oscillation 
times. Inside this radius stellar orbits react adiabatically 
to the forcing due to galaxy oscillations (see CD94), while 
outside galaxy oscillations become effective in changing 
significantly the stellar orbits. As expected, the effect be- 
comes however smaller and smaller for stellar orbits with 
increasing orbital times. Finally, note how there are not 
significant differences between the results in the large (Ca) 
and small (Cb) clusters: we interpret this as another man- 
ifestation of the comparable oscillatory periods in the two 
cases (see Fig. 1). 

6. Discussion and conclusions 

In this paper we presented a representative selection of 
numerical models aimed at the study of collisionless evap- 
oration of stars from cluster elliptical galaxies. From a 
theoretical point of view this possibility was pointed out 
in CG98, who showed that characteristic oscillation times 
of elliptical galaxies in near equilibrium configurations in 
the tidal field of the parent cluster are curiously compa- 
rable to the stellar orbital times in the outskirts of the 
galaxies themselves. From a numerical point of view the 
present work represents the natural follow-up of two recent 
preparatory works on the same subject (MC03ab), where 
we explored the simplest equilibrium configuration (i.e., a 
galaxy with the center of mass at rest at the center of a 
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models Cbi, Cbc, and Cbo, respectively. In all cases the maximum oscillation angle is set to tpM = 0.2 rad. 



triaxial cluster), but we left unaddressed the more compli- 
cate (although more astrophysically representative) case 
of galaxies with their center of mass in rotation around 
the (spherical) cluster center. We recall that in all the 
presented simulations the direct effect of the cluster gravi- 
tational field on stellar orbits was not included, in order to 
study the effects of galaxy oscillations only; we also recall 
that its effect is compressive inside the cluster core radius, 
while expanding outside (see, e.g., Valluri 1993, CD94). 
The main results can be summarized as follows: 

— A general analytical framework suitable for the study 
of orbital evolution of stars in galaxies oscillating 
around their stable equilibrium positions in the tidal 
field of the parent cluster is presented. In particular, 
we analyzed the cases of triaxial galaxies at the center 
of a triaxial cluster and in circular orbit around the 
center of a spherical cluster. 



— It is shown that, for realistic galaxy and cluster mod- 
els, galaxy oscillation times can be comparable to char- 
acteristic stellar orbital times in the external galactic 
regions, and thus important effects on stellar orbits are 
expected there. 

— Numerical simulations of galaxiesat the cluster cen- 
ter remarkably confirm the preliminary results pre- 
sented in MC03ab (where different galaxy models were 
adopted) . In particular, we found that orbits in the ex- 
ternal parts of the parent galaxy are affected by galaxy 
oscillations, with stars outside mi ~ 3 ~ 2.25i? nearly 
doubling the apocenters. Globally, a galaxy mass frac- 
tion of the order of 10% is affected by this process. 

— We speculate that the findings above could be related 
to the formation of cD galaxies. In fact, it is well known 
that numerical simulations of galactic cannibalism, al- 
though successful at reproducing BCGs, are in gen- 
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eral unable to obtain the extended halos of cDs. In 
this respect, the fact that the external parts of the 
galaxies are the only ones affected by collisionless evap- 
oration, is a nice property of the explored scenario. 
Unfortunately, our simulations cannot be used to study 
the resulting galactic profiles, and high-resolution N- 
body simulations should be performed to study further 
the matter. 

— For galaxies in circular orbits, the effect of collision- 
less evaporation is enhanced with respect to the pre- 
vious case, and stars can reach significantly larger dis- 
tances. This is not due to the direct effect of the CTF, 
but to galaxy oscillations only. In accordance with our 
preparatory analysis, we found that the radial trend 
of galactic oscillation periods at increasing distance 
from the cluster center is reflected by the location in 
the galaxy of the expanded region, being the most af- 
fected galaxies those orbiting near the cluster core ra- 
dius. Even in this cases, the affected galaxy mass frac- 
tion ranges is of the order of 5% - 10%, in interesting 
agreement with observational estimates. 

— The natural consequence of the results listed above is 
that, if collisionless evaporation is of importance in the 
production of the ISP (as direct contributor and/or 
help to the stripping phenomena due the direct CTF 
effect and to galaxy-galaxy encounters, as investigated 
in detail by Merritt [1983, 1984, 1985]), the density of 
ISP (when normalized to the local cluster luminosity in 
galaxies) should be enhanced around the cluster core 
radius, and rapidly declining at the cluster outskirts, 
due to the increase of the galaxy oscillation periods 
with the distance from the cluster center. 



Appendix A: Oscillation angles for a triaxial 
galaxy in circular orbit 

Here we give the explicit solution of eqs. (8) for assigned 
initial conditions ((po, ipo) and ($o, $o)- As shown in CG98, 
the two mixed equations of the second order 



-A-np - 
-A 2 §- 



B 2 <p, 



(A.I) 



(where the 4 coefficients A\ , A 2 , B\ , B 2 are given in eqs. 
[8]), can be separated in two identical biquadratic equa- 
tions, each of them involving only one of the two variables 
ip and For example, the equation for ip can be written 

as 



dip 
~dt T 

where 



A 



a tp 
dt 2 



B = 0, 



(A.2) 



A = A 1 + B X B 2 + A 2 



(u + v- l) 5 
uv 



+ ■ 



Q 2 C x 

, (4-3g)(l-n) 



B = A Y A 2 = fl 4 c x 



(4 -3«)(1 -«)(!-«) 



uv 



(A.3) 



(A.4) 



u = h/h, v = I 2 /h, and the functions f2 c (r) and 
q(r) are given in Sect. 3.2. The standard substitution 
ip = X exp(iut) in eq. (A.2) leads to the characteristic 
equation o> 4 — Auj 2 + B = 0, whose solutions can be writ- 
ten as 



A±; A± = 



A ± VA 2 - 4B 



(A.5) 



We finally conclude by pointing out that both orbital 
cases explored in this paper are rather exceptional. First 
of all, we only considered galaxy models that were ini- 
tially not tumbling nor rotating. In addition, most cluster 
galaxies neither rest in the cluster center nor move on 
circular orbits, but they move on elongated orbits with 
very different pericentric and apocentric distances from 
the cluster's center; in a triaxial cluster many orbits are 
boxes and some orbits can be chaotic. These latter cases 
can be properly investigated only by direct numerical sim- 
ulation of the stellar motions inside the galaxies, coupled 
with the numerical integration of the equations of the mo- 
tion of the galaxies themselves. In addition, a statistically 
significant galaxy population made by the sum of galax- 
ies with different masses, scale-lengths, and flattenings, 
should be considered. 
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note that in case of a stable equilibrium for the galaxy 
configuration (as assumed in this paper), the roots A± are 
positive, and that from eqs. (A.3)-(A.4) uj cx il c . Thus, 
the general solution of eqs. (8) can be written as 



p = acos(a>_t + a ) + bcos(ui + t + b ), 
= ccos(w_t + c ) + dcos(uj + t + do), 



(A.6) 



where the four amplitudes a, b, c, d and phases a , bo, Co, do 
must be determined by imposing the 4 initial conditions 
at t = 0. The remaining four constraints are provided by 
the request that eqs. (A.6) satisfy eqs. (A.l) at any time: 



c = a + it/2, 

do = b Q +7r/2; 

a{oj 2 _ - A)- cBilu- = 0, 

b{uj\ -A)- OIBXUJ+ = : 



(A.7) 



the identities involving amplitudes can be also written in 
terms of C and B 2 instead of A and B\. We can now 
proceed to impose the initial conditions on eqs. (A.6). The 
easiest way is to define the 6 new quantities 



r_ = a/c = Biw_/(u;2 - A), 
r + = b/d = Biuj + /{oj\ - A), 
A" = cos do, y = cos&o, 
, U = sinao, V = sin&o : 



(A.8) 
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note that r + and r_ are known. Evaluating eqs. (A. 6) at 
t = 0, and using eqs. (A.7)-(A.8) one obtains 

r-cX + r + dY = ipo, uu-cX + uo + dY = — #o, 

cU + dV = -i? , w_r_c[/ + u + r + dV = -^> , ( A - 9 ) 

X 2 + U 2 = l, Y 2 + V 2 = l. 

The first 4 equations above can be solved for the unknowns 
cX, dY, ell, dV in terms of r±, (d ,d ), (<p, ip ): then the 
two independent phases and amplitudes are obtained as 



c 2 = {cX) 2 + {cUf , d 2 = {dY) 2 + {dVf 
tanao = clI/cX, tan6o = dV/dY. 



(A.10) 



The remaining quantities are finally determined from eqs. 
(A.7)-(A.9). 

Appendix B: Potential, gravitational energy an 
TF components in homeoidal systems 

As in eqs. (11), (12) and (19) let use assume that 

p = p n p(m), (B.l) 

where m is defined as in eq. (13), and p n is a normalization 
constant for density. The total mass of the model (when 
finite) is given by 



M = 4iraia 2 a 3 p n M, 
where 



M 



p(m)m 2 dm. 



(B.2) 



(B.3) 



For example, in eqs. (11), (12), and (19) M = 1/(3 - 7), 
Ai = 1/30, and At = 7r/4, respectively. From potential 
theory, the general expression for <^ g (£) is found using the 
well-known formulas (Kellog 1953, Chendrasekhar 1969, 
BT87), 



r 

<X X ) = -Girp n aia 2 a 3 

Jo 

where 



A*[m(r, x) 
A(r) 



dT 



A(r) = ^/(a 2 + T)(a 2 +T)(a 2 +r), 



A*(m) = 2 / p(t)tdt, 

■>(r,x) 



and 



m: 



1=1 L 



(B.4) 

(B.5) 
(B.6) 

(B.7) 



For homeoidal density distributions the potential en- 
ergy tensor components are given by (Roberts, 1962) 



Uu = - 
where 



-Gtt 2 p^axa-iCLzalwi I A\P (m)dm, 
Jo 



uu 



aia 2 a 3 



dr 



o (a 2 +r)A(r)' 



(B.8) 



(B.9) 



and U = Uu + U22 + ^33- The three integrals Wi are 
calculated by using identities 238.03, 238.04 and 238.05 
from Byrd & Friedman (1971): 



Wi = 2a 2 a 3 



w 2 = 2d 2 a 3 



F(ip,k) - E(ip,k) 
k 2 sin 3 (p 

E(ip, k) - k l2 F{ip, k) - k 2 (a 3 /a 2 ) simp 



„. - (02/03) sin ip 
w 3 = 2a 2 a 3 — 



k 2 k' 2 sin 3 ip 
-E{fp,k) 



k' 2 sin ip 



(B.10) 



, (B.11) 



(B.12) 



where a 2 — a 2 /a 




a 3 /a\, ip = arcsin-^/l — a§, k 2 = 



(B.13) 



sin tdt, 



are the Elliptic Integrals of the first and second kind, re- 
spectively and k' 2 = 1 — k 2 . Also, we recall that the TF 
tensor is given by = —d 2 <fi/dxidxj-. the explicit evalu- 
ation of this quantity by using eq. (B.4) proves eq. (20). 
We also evaluate the quantity 2irGp c o needed there by us- 
ing the virial theorem, M c ay = —U. From eqs. (B.2) and 
(B.8) we obtain the exact expression 



27rGp n 



8Ma 2 



(B.14) 



For the density distribution in eq. (19) it can be proved 
that the value of the integral at denominator equals M = 
tt/4. 

Appendix C: Expansions for small flattenings 

The explicit form of the potential <f> generated by a generic 
homeoidal density distributions usually cannot be ex- 
pressed in explicit, closed form. On the contrary, its ex- 
pression in the limit of small flattenings is trivial. As de- 
scribed in Sect. 3, we introduce the two ellipticities with 
0-2/0-1 = 1 — e and a 3 /ai = 1 — r\\ obviously, in spherical 
models e = r\ = 0. We now expand eq. (B.4) for small 
flattenings, retaining only linear terms. We perform this 
expansion at fixed mass, i.e, we eliminate the coefficient 
iraia 2 a 3 p n between eqs. (B.3) and (B.4). With the nor- 
malization of lengths to 01, a tedious but straightforward 
calculation shows that <j> = —(GM/ai)4>, where 

M x0 = h{f) I (S + ^ + ^iM-Mr)] | 



f 3 
(e + r])h(r) _ (y 2 e + 5 2 i])h(r) 
3f 3 f 5 



(C.l) 



and I n (x) = f£ p(t)t n dt. Many dynamical properties of 
general density-potential pairs obtained with this trunca- 
tion procedure are presented elsewhere (for a full discus- 
sion of this technique see Ciotti & Bcrtin 2004). 
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h{f) 



For example, for ellipsoidal Hcrnquist models, 
f(2 + f) 



2(1 + f) 2 ' 



2(1 + ff 



T f(6 + 9r + 2f 2 ) o] . 



2(1 + f 



(C.2) 
(C.3) 
(C.4) 



Note that, at the center of the Hernquist model, the force 
remains finite but is not continuous, i.e., it depends on the 
approaching direction, and this fact produces undesidered 
orbital scattering due to amplification of numerical errors. 
This problem is avoided when using a density such that 
in eq. (12), from which 



h(f) 
h{f) 
h(f) 



P(10+10f + 5P+P) 
20(1 + f) 5 

f 3 (10 + 5f + r 2 ) 
30(1 + f) 5 ' 



5(1 + r) 



(C.5) 
(C.6) 
(C.7) 



Note that the force now vanishes at the model center. 

We also expand up to linear terms the three functions 
in cqs. (B.10)-(B.12) for fi -> and v -> 0: 



Wl ~ 2(1 - - i/) ( £ + ^ 



uto~2(l-/i)(l-i/)('± + ^±^ 



(C.8) 
(C.9) 
(CIO) 



With the aid of these expansions we can finally obtain 
also the linearization of the gravitational potential energy 
given in cq. (B.8) 



GM 2 J °° (m)rim / 
and eq. (B.14): 



+ 



2nGp, 



c0 



o?/™A$ (m)dm 



1 + 



/x + v 



2(n + u) 



(C.ll) 



(C.12) 
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